Traveling ion channel density waves affected by a conservation law 
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A model of mobile, charged ion channels embedded in a biomembrane is investigated. The ion 
channels fluctuate between an opened and a closed state according to a simple two-state reaction 
scheme whereas the total number of ion channels is a conserved quantity. Local transport mech- 
anisms suggest that the ion channel densities are governed by electrodiffusion-like equations that 
have to be supplemented by a cable-type equation describing the dynamics of the transmembrane 
voltage. It is shown that the homogeneous distribution of ion channels may become unstable to 
either a stationary or an oscillatory instability. The nonlinear behavior immediately above thresh- 
old of an oscillatory bifurcation occuring at finite wave number is analyzed in terms of amplitude 
equations. Due to the conservation law imposed on ion channels large-scale modes couple to the 
finite wave number instability and have thus to be included in the asymptotic analysis near onset of 
pattern formation. A modified Ginzburg-Landau equation extended by long-wavelength stationary 
excitations is established and it is highlighted how the global conservation law affects the stability 
of traveling ion channel density waves. 

PACS numbers: 87.10.+e, 05.65.+b, 87.16.Uv, 47.20.Ky 



INTRODUCTION 



Spatiotemporal pattern formation is omnipresent in 
physical, chemical as well as biological systems driven out 
of equilibrium and can be be ascribed to a finite number 
of universality classes 0, 0, 0, ■ Cardiac spiral 
waves as well as the propagation of nerve signals along 
the axon of a neuron are prominent exam ples of bioelec- 
tric nonstationary patterns |E 0, 0, ITU Il2l liljj : other 
examples might include the emer gen ce of aster-like pat- 
terns in filament-motor-systems 0, ^E 03 or the 
nematic ordering of biopolymers such as actin and mi- 
crotubule complexes Hij- Futhermore spontaneous 
pattern formation of ion channels is well-known to oc- 
cur in cell membranes ,2Gj, 2 1 . 22 : 23] that are a ubiqui- 
tous building block in biology. In most biological systems 
the pattern forming processes are more elaborate than in 
fluid dynamical ones. Hitherto a qualitative, especially 
a quantitative understanding of universal aspects of pat- 
tern formation has mainly been achieved in fluid systems 
[E0]< Biological systems, such as the model investigated 
in the present work, give however rise to new universal 
features of pattern formation. The global conservation 
law imposed upon the number of ion channels embedded 
into the cell membrane leads within the considered model 
system to a modified set of amplitude equations describ- 
ing the generic properties of traveling ion channel density 
waves and alters moreover their stability behavior. 

The considerable amount of information on membranes 
gained throughout the years has been synthesized into 
the fluid-mosaic model of a biomembrane |24| which in- 
cludes the following main features: the phospholipids 
serve as a solvent for proteins and as a permeability reg- 
ulator, the latter by adopting an energetically very ef- 
fective bilayer configuration; membrane proteins are free 
to migrate within the lipid-bilayer; the conductance of 
the membrane layer is mainly made up of discrete ion 



channels formed by protein molecules [25| . Freely mov- 
able ion channels embedded in a fluid-mosaic membrane 
have been observed to self-organize if the concentration 
of salt across the membrane exceeds a certain threshold 
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The basic fluid-mosaic model underwent however con- 
siderable refinement as it neglects first of all that ion 
channels may fluctuate between opened and closed states 
due to interactions with signal molecules [2f| and sec- 
ondly that channel proteins may become immobilized 
due to binding to the cytoskeleton |3(|. A compound 
system of a membrane with embedded ion channels is 
typically driven out of equilibrium due to ionic concen- 
tration gradients and transmembrane fluxes: Opened ion 
channels induce ionic transmembrane currents looped to 
the transmembrane voltage. Therefore thermodynamic 
models have to be extended by ion transport mechanisms 
based upon Nernst-Plank theory and by electrodynami- 
cal aspects |2E |2^ ■ The resulting electrodiffusive models 
can be realized by considering a membrane seperating a 
narrow cleft of electrolyte from an electrolytic bath of 
ions in, a configuration often met in biological situa- 
tions: closely neighboring cells like post-synaptic mem- 
branes of a neuronal synapse or membrane cables forming 
dendrites and axons of neurons are important examples 
of such geometries. The transmembrane voltage corre- 
sponding to the voltage difference across the thin layer 
of electrolyte is then given by a one- or two-dimensional 
cable- type equation |2lt l32| . 

It has been shown that the coupled dynamics of chan- 
nel proteins and a binding-release reaction of the latter 
ones with the cytoskeleton leads to spatial peridoc pat- 
terns m described near onset of pattern formation by a 
Ginzburg-Landau equation [jj^j . A model of two differ- 
ent, noninteracting species of ion channels has been stud- 
ied in Refs. [3411351 ] A generic evolution equation for non- 
stationary patterns generated by two different ion chan- 
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nel species assumed to remain permanently in an opened 
state was given in Ref. 35] wherein a constant number 
of opened ion channels was presumed. In the present 
work we will consider ion channels that may switch be- 
tween an opened and a closed state 1361, |37j according to 
a simple two-state reaction scheme with constant rates, 
the diffusion coefficients as well as the effective charges of 
the channel proteins being state-dependant. If compared 
to the models discussed in Refs. |35|, the system at 
hand may be interpreted as one comprising two types of 
ion channels with a reaction stimulating transition be- 
tween the different channel species while imposing the 
constraint of a conserved total number of ion channels. 

After introducing the model in Sec. [H] the linear insta- 
bility analysis performed in Sec. IIIII evinces that either 
a stationary or an oscillatory bifurcation is met and is 
concluded by a phase diagram in Sec. lIIIDl Immediately 
above the threshold of a supercritical bifurcation the am- 
plitudes of the unstable pattern forming modes emerging 
from the homogeneous basic state are small enough to ap- 
ply the powerful perturbational technique of amplitude 
equations. Oscillatory instabilities have been broadly dis- 
cussed in literature as long as systems lack global conser- 
vation laws and their nonlinear evolution is governed by 
the most prominent class of amplitude equations, the so- 
called complex Ginzburg-Landau equations 0, Ei] ■ The 
influence of a conserved quantity on pattern forming in- 
stabilities, whether stationary or oscillatory, has been, 
starting from symmetry arguments, a subject of continu- 
ous interest during the recent years [3^.l38ll3^.l40ll4lll42] | 
but was never discussed for a specific microscopic model. 
As stationary bifurcations have already been studied for 
the present biological system in Ref. |43j |. we will focus 
in the remaining sections on Hopf-bifurcations to travel- 
ing ion channel density waves. Within the context of a 
weakly nonlinear approach near onset of pattern forma- 
tion in Sec. II VI the complex Ginzburg-Landau equation 
known to describe stripe patterns is altered through the 
coupling to an evolution equation for a large-scale, real 
mode evoked by the presence of a global conservation 
law as shown in Sec. IIV Al In Sec. IIV CI we investigate 
how the implication of the large-scale mode reflecting the 
conservation law modifies the stability regions of super- 
critically bifurcating traveling ion channel density waves 
obtained in Sec. IIV Bl It turns out that the latter ones 
become always unstable with respect to an instability at 
finite wave number and that observing stable stripe pat- 
terns becomes less probable the more ion channels con- 
tribute to the total conductance of the membrane layer. 
Analytical calculations are confirmed by numerical simu- 
lations in Sec. IIV Dl Concluding remarks and an outlook 
on future work are given in Sec. 



II. 



MODEL SYSTEM 



We consider a model system constituted of a membrane 
and embedded ion channels which may be encountered 



either in an opened or in a closed state. The dynami- 
cal equations of opened (o) and closed (c) ion channel 
densities n (r, t) respectively n c (r,t) are coupled via a 
global conservation law J dS{n + n c ) = const, (wherein 
dS denotes the area of the considered biomembrane) im- 
plemented by a simple local reaction scheme as well as to 
the cable equation governing the evolution of the trans- 
membrane voltage. 



A. Basic equations 

The switch between the opened and closed state of 
an ion channel, referred to as opening-closing reaction, 
is governed by the following monomolecular reaction 
scheme 



n c ^ n Q 



(1) 



with constant rates r c and r . According to the homoge- 
neous, stationary ion channel densities 
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(2) 



channel dynamics are best captured in terms of the de- 
viations n n and n r 



n Q n n , 



n r = n r . - n c 



(3) 



from their respective homogeneous mean densities n D , n c . 

Local transport mechanisms as well as the mentioned 
conversion kinetics suggest the electrodiffusion-like equa- 
tions of motion to be 

d t n = -Vj - r c h + r h c , (4a) 
d t n c = -Vj c -I- r c n a - r D n c , (4b) 

wherein the current densities j DiC depend on the chemical 
potential /i 0jC and the electrophoretic drift |20| via the 
transmembrane voltage v as follows 



- V/x 0iC - v ^ c q , c n , c Vv . 



(5) 



The occuring mobilities v D ,c of opened and closed ion 
channels 



D„ 



(6) 



are proportional to the respective diffusion constants 
Z? 0iC . The stated charges q 0tC of opened and closed ion 
channels are considered to be effective charges taking 
screening as well as electroosmotic effects into account 
|22j |. If excluded volume interactions between opened 
and closed channels together with nonlinear effects are 
neglected, the chemical potential takes the simple form 



(7) 
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Presuming a one-dimensional system substitution of the 
current j o c according to Eq. (J5J in Eqs. 10} gives rise to 
the following equations 

d t h = D \7 2 n + ^° q ° X7 x \{n + n a )\7 x v] - r c h Q + r h c 

(8a) 

r c n 
(8b) 



d t h c = D c \7 2 h c + -r-^V ' x [(h c + n c )\7 x v] - r a n c + r c h. 



describing the spatiotemporal evolution of the ion chan- 
nel densities. Those equations have nevertheless to be 
supplemented by an equation governing the transmem- 
brane voltage v. 

As a biomembrane is electrically equivalent to a cable 
whose repetitive basic unit is formed of a capacitor mod- 
eling the pecularities of lipid molecules, a voltage source 
reflecting the properties of an opened ion channel and 
a leak conductivity taking all remaining effects into ac- 
count, the Kelvin cable equation may be considered to 
determine the transmembrane potential v on condition 
that the caracteristic length scales of lateral patterns are 
large compared to the width of the electrolyte layer. A 
derivation of the one-dimensional cable equation 



C m d t v 



1 



^V 2 u 



A n (v - E a ) - Gv . 



(8c) 



from the general Nernst-Plank theory was given in |32l| . 
Herein C m denotes the membrane capacitance, fi the re- 
sistance of the thin membrane layer, G the leak conduc- 
tance and A n the conductance of the opened ion chan- 
nels per unit length. Finally E Q stands for the Nernst po- 
tential, i.e. the resting potential of the membrane layer. 

In most cases we consider the limit of a very quick 
response of the transmembrane voltage, i.e. formally 
C m = holds, which is characteristic for physiological 
conditions 361. 



B. Scaled equations 

We rescale Eqs. (|SJ) by introducing dimensionless co- 
ordinates for space x' = xj A and time t' = t/r with 
a typical length scale of an electrical perturbation A = 
[Q(A n + G)]" 1 / 2 and the time constant of displacement 

2 

t = A / D a . As normalized particle densities we opt for 



Tlr , , 

=r> 9 

Tin 



the rescaled voltage reading V :— (v — vn)q /kBT, 
wherein vr — aE Q denotes the resting voltage and 
a = A n /(A n + G) € [0,1] the density parameter 
measuring to what extent the ion channels contribute to 
the total conductance of the membrane layer. Through 
defining a normalized relaxation time ry = ttC m D , we 



obtain on the basis of the additional abbreviations 

q E 



e = — 



R = 



<lo 



A 2 


D = 


D c 








0O = TT , 


f3c = 


fr c 



(10b) 



the following scaled equations 

d t N Q = ^JV + V X [(1 + JV )V X V] 

-f3 c (N - N c ) , (11a) 
d t N c = DV 2 x N c + DRV x [{l + N c )VV] 

+/3 (N -N C ), (lib) 
T V d t V = (V£ - 1)V - a(l - a)eN - aiV V(llc) 

wherein the primes of the new coordinates x' and t' have 
been suppressed for the sake of simplicity. 

It should be noticed that the spatially homogeneous 
basic state of these scaled equations is N a = N c = V = 0. 



III. INSTABILITY OF THE HOMOGENEOUS 
STATE 

In a certain parameter range the spatially homoge- 
neous basic state, namely N Q = N c = V — 0, becomes 
unstable against infinitesimal perturbations and a spa- 
tially inhomogeneous pattern is likely to emerge from this 
very state. In order to determine this interval of instabil- 
ity we transform the linear part of the partial differential 
equations Eqs. i|llfl by means of the ansatz 



N \ 








= ( £ j e At+2k r 


(12 


v ) 







into a set of homogeneous linear equations for the am- 
plitudes N , N c and V. These equations can be fur- 
ther simplified by assuming an instantaneously adapt- 
ing transmembrane voltage V, i.e. presuming ry = 

{C m w 1^, ^ » 10 8 Ohm, D w 0.1^ HI) or equiv- 
alently C m — 0. Within this limit the voltage perturba- 
tion V may be expressed in the linear regime in terms of 
the opened ion channel density N a 

with e playing the role of the control parameter 

e = a(l-a)e, (14) 

and may thus be adiabatically eliminated from the other 
equations. For the two remaining linear and homoge- 
neous equations governing the pure ion channel dynam- 
ics 



-Pa 



DRek 2 



= 0, 



(15) 
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the solubility condition provides a quadradic polynomial 
for A 



X + AX + k B = . 



the two coefficients A and B reading 



A = 



p Q + f3 c + (1 + D)k 2 - 



ek 2 



B = f3 + D(f3 c + k 2 ) - 



1 + k 2 ^ 

0o + D{(3 C R + k 2 ) 
1 + k 2 



(16) 



(17a) 



(17b) 



Depending on the model parameters the two solutions 
Ai,2 of Eq. (|16|) are either real or complex. Provided that 
the real part of only one of those solutions Ai,2 becomes 
positive in the neighborhood of a wave number k, the 
basic state N = N c = V — is said to be unstable: In 
case of a vanishing imaginary part Im(A) one has to deal 
with a stationary instability 2] whereas for finite imagi- 
nary parts, i.e. u> = Im(A) ^ 0, an oscillatory instability, 
a so-called Hop] "-bifurcation Q, is met. 

Forthcoming discussions are considerably clarified by 
introducing the conversion rate ratio 

/?=§ = -• (18) 
Pc r c 



A. Dispersion relations 

The real parts o\ t i = Re(Ai,2) of Ai,2(fc), sketched in 
Fig. ^for D — 1 and (3 = 0.4, exhibit typical shapes as a 
function of the wave number k. Beyond the threshold of 
pattern formation these different shapes of o~i^{k) may 
lead in the weakly nonlinear regime to different univer- 
sality classes of amplitude equations as discussed in the 
following paragraphs. 

It should be noted that for typical Nernst potentials 
with absolute values reaching from \E Q \ = to \E„\ = 100 
mV (equilibrium potentials of the most important volt- 
age gated ion channels for the mammalian skeletal mus- 
cle: E a = —98 mV in case of K + channels and E Q = +67 
mV in case of Na + channels (3(J), the control parameter 
e defined in Eq. i|14|) takes values in the range |e| G [0, 10] 
under physilogical conditions. Thus all instabilities pre- 
sented in Fig.^may play a role in a real biological system 
as electrophoretic charges are usually of a few elementary 
units. 

Since the number of ion channels is a conserved quan- 
tity for the investigated model, i.e. J dx(N Q + N c ) = 
const., the largest real part labelled eri(fc) is diffusive in 
the long-wavelength limit 



Ai(fc -> 0) = <7i(fc -> 0) 



-D N k 2 



0(fc 4 ) 



with D N = [Dp c (Re - 1) + (3 a (e - 1)] / (0 O 
confirmed by all examples shown in Fig. ^ 



(19) 

Pc) as 



In part (a) of Fig. the real parts of both eigenvalues, 
(71,2 (fc), are drawn as a function of the wave number k 
for a parameter combination that leads to a stable homo- 
geneous state. Part (b) shows an example for the same 
charge and conversion rate ratio R = 1.1 respectively 
P = Pol P c = 0.4, with a positive curvature of the real 
part in reference to the wave number k 

Xi(k -> 0) = ai{k -> 0) = D N k 2 - 0.9fc 2 > (20) 

and a maximum of the dispersion relation at a finite 
value k m ~ 0.6. If the charge ratio R is reduced while 
the conversion rate ratio /3 remains the same, the homo- 
geneous state may become stable with respect to long- 
wavelength perturbations as for the parameter set used 
in part (c), whereas the curvature <9/c<7i(fc)| fe=0 < is 
negative. The homogeneous state gets nevertheless un- 
stable against perturbations with a finite wave number 
in the vicinity of k = k c = 0.8 as can be foreseen from 
the extremum condition dkO~i(k)\ k=k = 0. 

As the charge ratio R is further diminished, the two 
real parts ci^k) approach each other as depicted in part 
(d) . Decreasing R even further leads to a pair of complex 
eigenvalues Ai ; 2 instead of real ones [cf. Fig.[T](e)], mean- 
ing that the homogeneous basic state becomes simulta- 
neously unstable against a pair of traveling waves (solid 
branch) with a critical wave number of k ~ 0.8 and a spa- 
tially periodic modulation (dashed branch) with a wave 
number of k ~ 1.35. As R gets even smaller the latter 
stationary branch, i.e. the dashed one visible in part (e), 
becomes strongly negative and only a Hopf-bifurcation 
to traveling waves remains as depicted in part (/). When 
one of the dispersion curves in Fig.Q]hits the k— axis as a 
function of the control parameter e, the threshold of the 
respective instability is met. 



B. Stationary threshold 

The threshold ef for a stationary bifurcation follows 
from the neutral stability condition a = A = 0, which is 
equivalent to claiming B = in Eq. I|16|l . From B = 
the wave number dependent neutral curve Q, 

_ [D (k 2 + p c ) + p o ] (1+g) 
eW ~ D{k 2 + R p c )+p o {21 > 

may be derived. This neutral curve takes its minimum 
either at a vanishing critical wave number = 



e S (0) = 



DPc + Po 

DRP C + Po 



(22) 



or the minimum becomes even lower than the latter one 
for a finite value k^ of the wave number 



/?o , \/Pc{R — 1) [D (RPc — I) + Po] 
(k c ) --HPc -+ -j= , 

(23) 
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FIG. 1: The real parts of the two solutions Ai,2(&) of Eq. I16H 
are shown as a function of the wave number k for different 
values of the charge ratio R while the conversion rate ratio 
13 = Po/Pc = 0.4 with f3c = 0.5 and the diffusion ratio D = 1 
are kept constant. Dashed branches correspond to stationary 
bifurcations, while solid ones to oscillatory instabilities. In 
part (a) the charge ratio is R = 1.1 and the control param- 
eter has a subcritical value e s = 0.93, whereas e is always 
supercritical in the remaining cases. In part (&) the ratio 



R 



1.1 and e s = 1.77, in (c) R = -0.1 and e s 



3.04, in 



part (d) R = -0.39 and e b = 3.69, in part (e) R = -1.244 and 
e H = 5.32 and finally in part (/) R = -1.8 with e H = 5.83. 

where kf is given by the extremum condition dte (k) = 
0. The threshold value of the control parameter e at kf 
reads 

s = [DPc{R-1)-T] [D(Rp c -l)+po-T] 

DT 1 ' 

wherein the abbreviation 

T = y/D/3 c (R - 1) [D {Rf3 c - 1) + f3 ] (25) 
has been introduced. 

C. Hopf-bifurcation 



0.0 0.6 1.2 1.8 u 0.0 0.2 0.4 0.6 0.8 1.0 k 

FIG. 2: Neutral curves for a stationary (e s ) and oscillatory 
instability (e H ) as a function of the wave number k (on the 
left hand side) and the frequency along these neutral curves 
(on the right hand side) as given by Eq. 12711 for the parameter 
set used in Fig.^e). 



in order to satisfy the solubility condition Eq. I|lbf) . they 
provide two conditions allowing to determine the neutral 
curve of the oscillatory instability 

eg(fc)= (l + * a )(* a +^+ft> + ft) > (26) 

as well as the frequency 

oj 2 (k) = k 2 B(k,e) (27) 

with B from Eq. (|17b|) . Both quantities are shown as 
a function of the wave number k for a typical parame- 
ter combination in Fig. [21 where the stationary and the 
Hopf-branch have a similar threshold as expected from 
Fig. de). If parameters corresponding to Fig. WLf) are 
chosen, the shape of the neutral curves depicted in Fig. [2 
would not be altered but the Hopf-branch would have the 
lowest threshold. The neutral curve specified in Eq. I|26() 
takes its minimum at 



k H = 



1 + D 



(28) 



and the corresponding critical value of the control pa- 
rameter at k? is 



whereas the Hopf-frequency evaluates to 



(29) 



k^J-fVTTD - Dk" 2 (D + 20 o ) (30) 



with 



7 = PoVPo + Pc + RPce? 



(31) 



At onset of a Hopf-bifurcation the real parts 01,2 of the 
complex conjugate pair of eigenvalues Ai^ vanish, while 
the imaginary parts lo\^ = ±Im(A) remain nevertheless 
finite. Assuming complex eigenvalues A1.2 the polynomial 
in Eq. (|l(j|) may be decomposed into its real and imag- 
inary part. Since both parts have to vanish separately 



D. Neutral curves and phase diagram 

There are obviously three different types of instabili- 
ties of the homogeneous state occuring in different areas 
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FIG. 3: The figure shows the respective ranges of a stationary 
bifurcation either of the Cahn-Hilliard-type with a critical 
wave number k c = 0, uj c — or with a finite critical wave 
number k c ^ 0, lu c = and of a Hopf-bifurcation fc c 7^ 0, uj c 7^ 
0. The symbols (a) to (/) correspond to the parameters in 
the i?-/3-plane for whom the dispersion relations are plotted 
in Fig.Q parameters being D = 1 and j3 c = 0.5. 

of the i?-/3-plane, /? being the conversion rate ratio de- 
fined in Eq. I|18|). One has either a Cahn-Hilliard-likc 
|44| instability characterized by e s (k = 0) < ef and 
e s (fc = 0) < e^ (III s -bifurcation 0), i.e. the long- 
wavelength excitations dominate; a stationary instability 
at a finite wave number k with ef < e s (k = 0) as well as 
ef < e^ (I s -bifurcation [2]); or a Hopf-bifurcation distin- 
guished by < e s (k — 0) and < ef (I Q -bifurcation 
2] ) . These three regions of existence are shown in Fig. [3] 
for (3 C — 0.5 and D = 1, wherein the dashed line is deter- 
mined by the condition that the two stationary thresh- 
olds agree 

ef = e s (fc = 0), (32) 

whereas the solid line is given by the codimension-2 con- 
dition 

ef = ef. (33) 

The dispersion relation corresponding to the point la- 
belled (c) in Fig. |3 is given by Fig.^c). The parameter 
set corresponding to (e) in Fig. [3] lies in the immediate 
vicinity of the codimension-2 line and the neutral curves 
as well as the change of frequency along those curves are 
drawn as a function of the wave number k in Fig. [21 the 
associated dispersion relation being depicted in Fig.^e). 

If the diffusion ratio D is reduced compared to D = 1 
in Fig. 13 while the conversion rate (3 C is nevertheless kept 
constant at /3 C = 0.5, the area wherein a I s stationary bi- 
furcation occurs is mostly affected as depicted in Fig. 0] 
For large values of the diffusion ratio, such as D — 0.8 in 
Fig-EI a )j a stationary bifurcation at finite wave number 




1 2 3 4 a 5 1 2 3„4 



FIG. 4: The boundary between a I s -bifurcation and a Ills- 
bifurcation (dashed line) as well as the one corresponding to 
the transition from a stationary to an oscillatory bifurcation 
(solid line) are shown for different diffusion ratios D — 0.8 in 
(a) and D — 0.18 in (b) while the conversion rate /3 C = 0.5. 



is encountered in a broad band of the i?-/3-plane bounded 
by the dashed line corresponding to the border between 
the III S - and I s -bifurcation and the solid line visualizing 
the codimension-2 line. Reducing D induces the region of 
existence of the I s -bifurcation to shrink significantly. For 
very small diffusion ratios, such as D = 0.18 in Fig.^Jb), 
a broadening of the I s -bifurcation-region occurs for large 
values of f3 or equivalently for strongly negative values of 
R, which is not so pronounced for intermediate diffusion 
ratios D, while its area of existence remains strongly con- 
fined for slightly negative charge ratios. For 1 < j3 < 3 
the solid line seperates in case of D = 0.18 [Fig. ^b)] a 
Hopf-bifurcation from a III s -bifurcation instead of a I s - 
bifurcation otherwise, such an interval of /3 values can 
always been determined as long as D < 0.3. 

Furthermore it should be noticed from Fig.0|that with 
decreasing values of D the slopes of the codimension- 
2 line and of the phase limit between the III S - and I s - 
bifurcation become steeper which progressively shifts the 
domain of emergence of an oscillatory bifurcation, lying 
beneath the codimension-2 line, to more negative charge 
ratios R. 

If the diffusion ratio D is however kept constant while 
the conversion rate (3 C is varied, the relative positions 
of the phase boundaries remain almost unchanged as it 
may be observed in Fig. [S] wherein D = 1. The latter 
figure evinces that the area of occurence of a I s stationary 
bifurcation located in-between the two phase boundaries 
gets insignificantly smaller for increasing values of (3 C . 
Hence the value given to the conversion rate (3 C is not 
determinant while focusing on Hopf-bifurcations in the 
forthcoming sections. 



IV. WEAKLY NONLINEAR ANALYSIS 
BEYOND A HOPF-BIFURCATION 

We will focus in the present section on the oscilla- 
tory bifurcation regime while a partial analysis of the 
I s -bifurcation may be found in Ref. 01 • Immediately be- 
yond the threshold of a supercritical bifurcation, where 
the amplitude of a pattern emerging from the homoge- 
neous state is still small, the concept of the so-called am- 



7 



A. Generic amplitude equations 




FIG. 5: Phase diagram in the i?-/3-plane for different values of 
the conversion rate /3 C given by Eq. 1181 . The codimension-2 
lines (solid lines) as well as the phase boundaries between a I s - 
bifurcation and a III a -bifurcation (dashed lines) are displayed 
for (3 C = 0.3, /3 C = 0.5 and f3 c = 0.8. The directed arrows 
visualize the shrinking of the area wherein a I s -bifurcation is 
encountered with growing values of the conversion rate /3 C . 
The diffusion ratio D is assumed to evaluate to 1. 



plitude equations is a very successful one to characterize 
the nonlinear behavior of patterns as it is exemplified for 
several physical, chemical and biological systems in Refs. 

Due to the presence of the global conservation law, 
namely J dx(N a + N c ) = const., all dispersion re- 
lations depicted in Fig. ^ behave diffusively in the 
long- wavelength limit. Therefore large-scale modes are 
marginal at the considered bifurcation point and have 
to be included in the reduced dynamics near onset of 
pattern formation. Restricting our investigation to trav- 
eling ion channel density waves, the complex Ginzburg- 
Landau equation |3^| normally used to describe travel- 
ing waves has to be extended by an evolution equation 
for the large-scale mode that is evoked by the aforemen- 
tioned global conservation law. This set of coupled am- 
plitude equations, is obtained in Sec. IIV Al by a system- 
atic perturbation expansion from the scaled equations 
Eqs. whereby the coefficients occuring in these am- 
plitude equations reflect the dependance of the pattern 
on the parameters of the underlying system. The sym- 
metry consitstency of the latter equations is checked in 
Sec. HVAl as well. 

After having studied the bifurcation structure by 
means of the complex Ginzburg-Landau equation in 
Sec. IIV Bl it will be shown in Sec. IIV 01 that the implica- 
tion of the large-scale mode dramatically changes, in case 
of an infinitely extended system, the results on stability 
of traveling waves acquired in Sec. IIVBI Finally the va- 
lidity of the performed amplitude expansion is checked 
numerically in Sec. IIV Dl 



The small parameter used in the perturbative deriva- 
tion of the coupled amplitude equations is the relative 
distance rj to the threshold e c or likewise e c 



V 



(34) 



Given a system presenting a global conservation law, 
as encountered in this work, the generating field wi = 
(N ,N C ) T of traveling waves in the vicinity of the bifurca- 
tion point comprises a complex, propagating mode A and 
a real, large-scale mode C stipulated by the constraint of 
conserved number of ion channels 



Wl = Ae e i( ^ t+k "- r) + C*v 



+ c.c. 



(35) 



wherein c.c. denotes the complex conjugate of the pre- 
ceding terms. Comparable to the adiabatic elimination 
of the transmembrane voltage V performed via Eq. i|13|) 
in Sec. II I II within the linear regime, V can be similarly 
expressed in terms of the channel densities within the 
weakly nonlinear regime, the field wi becoming thus im- 
plicitly dependant of the membrane voltage, in the 
latter ansatz designates the critical wave number defined 
in Eq. (|28[1 and the cited eigenvectors are 



e 



1 

F 



vo 



(36) 



the second component of the eigenvector eg reading 



Fa 



(37) 



The amplitudes A and C of the fluctuations Wi depend 
on the following critically slow time scales T\ — y/rjt and 
T = nt as well as on the scaled spatial variable X — 
^Jrjx, thus suggesting a standard multiscale perturbative 
approach |2| , that has to be combined with the aforesaid 
adiabatic elimination of the transmembrane voltage as 
sketched for convinience in Appendix lAl 

The expansion of the fundamental equations Eqs. Hll|) 
and of the pattern forming fields N , N c as well as V 
with respect to powers of the small parameter ^/r j results 
on the basis of the ansatz specified in Eq. (|35[1 in the 
subsequent system of coupled amplitude equations 

t (d t + v g d x ) A = (77 + io) A + £g (1 + ic ) d 2 x A 

~ gi (l + i Cl )A\A\ 2 -g 2 {l + ic 2 )AC 
- 53 (1 + ic 3 ) C 2 A + 04 (1 + ic 4 ) Cd x A 
+ g 5 (1 + ic 5 ) Ad x C , (38a) 
d t C = D x dlC-bd x \A\ 2 , (38b) 

provided that space and time are once more rescaled to 
x and t. 

The analytical expressions of the coefficients appear- 
ing in the latter equations Eqs. (|38|) as functions of the 
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parameters of the underlying biological model are rather 
lengthy. Instead of giving their analytical forms we there- 
fore plot them in Fig.EJas functions of the conversion rate 
ratio (3 while keeping the remaining parameters constant. 

Likewise amplitude equations, wherein the real, how- 
ever damped, mode C plays the role of a concentra- 
tion field, were used to describe localized traveling-wave 
trains observed in experiments on binary mixture con- 
vection and to clarify the influence of a long-wave 
mode on solitary waves [46| . The amplitude equations 
discussed in the latter references were recovered by a 
reaction-diffusion system |4^. In Ref. the nonlin- 
ear behavior of traveling waves arising in a supercriti- 
cal Hopf-bifurcation, coupled to a real, weakly damped 
mode C that is advected by the waves and thus affects 
the stability of the latter ones, yielded envelope equations 
resembling to Eqs. (J2HJ- I n a h the previous examples the 
real mode C was nevertheless damped and not a true 
zero mode, i.e. a conserved quantity: Slightly different 
amplitude equations, motivated by work on magnetocon- 
vection as well as rotating convection, with a true con- 
served quantity were derived by symmetry arguments in 
Ref. ji^l . The amplitude equations of Ref. [43 as well 
as those found for the present biological model will be 
further discussed within the scope of a two-component 
reaction-diffusion model which also describes the cou- 
pling of a traveling wave to an oscillatory long-wavelength 
mode instead of a stationary one and leads thus to a gen- 
eralized version of Eqs. I|38[l |49| . 

Let us now check if the amplitude equations Eqs. i|38|) 
obtained by the multiscale perturbative approach are 
consistent with the symmetries of the underlying model 
system. The electrodiffusive model equations (JSJ or 
equivalently the scaled basic equations {TJ| present trans- 
lational as well as reflectional symmetry. Accordingly 
amplitude equations describing standing waves generated 
by two complex counterpropagating linear modes A\ (x, t) 
and A 2 (x,t) coupled to a real, large-scale field C{x,t) 
evoked by the conservation law inherit the following in- 
variance under spatial and temporal translation: 



x^x + S: A l -> A l e l6 1 A 2 A 2 e-' l \ C -> C, 

(39a) 

t^t + n: A 1 -► A 1 e ict >,A 2 -> A 2 e^, C^C, 

(39b) 



where 9 = k^S, d> = u r n as well as under reflection: 



x 1— > —x : Ai — > A 2l A 2 — ► At, C -> C . (40) 



Due to the scalar constraint of a total conserved num- 
ber of ion channels the translational symmetry, c.f. 
Eq. I|39a|) , stipulates that the leading-order coupling term 
that extends the complex Ginzburg-Landau equations 




FIG. 6: The coefficients of the coupled amplitude equations 
Eqs. 11381 are plotted as functions of the conversion rate ratio 
j3 for the parameter combination D — 0.5, a = 0.8, R = —5.5 
and p c = 0.5. 



r(d t ±v g d x ) Ai 



(rj + ia) Ai + $(l + ico)8£Ai 
- g 1 (l + ici)A i \A i \ 2 
-geil + ic^A^AA 2 , 

with i,j G {1,2} andj ^ i , (41) 

is proportional to Ai C (i=l,2). Possible higher order 
terms that might be included in extended Ginzburg- 
Landau equations comprise AiC 2 or C(d x Ai) + Aid x C 
with i = 1, 2. In order to satisfy the reflectional symme- 
try given by Eq. (|40|l the complex coefflcents of the latter 
two terms need to be antisymmetric with respect to re- 
flections as they involve first order spatial derivatives. 

In absence of eventual coupling terms the real mode 
C fulfills a diffusion equation as suggested by the con- 
served quantity and the reflectional symmetry speci- 
fied in Eq. (|40|l . The translational symmetries cited in 
Eqs. H39|) then claim that coupling terms from the prop- 
agating modes A± and A 2 likely to occur in the evolution 
equation of C involve equal numbers of Aj and their as- 
sociated complex conjugates A* (£ = 1,2). Therefore the 
leading coupling term consistent with reflectional symme- 
try reads c^d^il 2 — l^l 2 ), the generic evolution equation 
for C finally being a diffusion equation with nonlinear 
forcing from Aj (i = 1, 2). 

Within the limit A\ — > A and A 2 — > 0, the standing 
wave patterns reduce to traveling waves and the symme- 
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try consistent amplitude equations happen to be exactly 
those, namley Eqs. i|55)l. obtained by a standard multi- 
scale approach from the inital model equations Eqs. 111(1 
in the previous section. Furthermore the coefficients 54, 
55 and b appearing in Eqs. l{3*5ll are pseudoscalars as re- 
quired by the reflectional symmetry of the underlying 
model system. 



B. Bifurcation structure of traveling waves in 
absence of the large-scale mode C 

Since the large scale mode C only affects stability prop- 
erties of traveling waves, we may neglect C while deter- 
mining their general bifurcation structure. The problem 
at hand then corresponds to a simple traveling wave of 
complex amplitude A and group velocity v g described by 
the amplitude equation 



r (d t + v g d x ) A = (rj + ia) A + f§ (1 + ic ) 
- gi (1 + ia) A\A\ 2 . 



diA 



(42) 



By changing to a coordinate frame moving with the 
group velocity v g , i.e. by performing a Gallilei trans- 
form, Eq. H42|) can be reduced to the prominent complex 
Ginzburg-Landau equation 0, 

d t A = (77 + ia) A + (1 + ic ) d 2 x A -gi(l + id) A\A\ 2 , 

(43) 

wherein the linear imaginary part ia is easily eliminated 
by the transformation A — > Ae lat . 



1. Supercritical vs. subcritical bifurcation 

For the complex Ginzburg-Landau equation it is well 
known that traveling waves bifurcate supercritically if 
the real part of the coefficient of the cubic term A\ A\ 2 in 
Eq. (|43[1 is positive, namely gi > 0. The tricritical line, 
along which the coefficient g\ vanishes, marks the bor- 
der between super- and subcritically bifurcating waves 
in parameter space 33] . 

Fig.0shows how the tricritical point, namely the root 
of the coefficient g±, moves at a strongly negative charge 
ratio R — —5.5 towards greater values of the conversion 
rate ratio P as D becomes larger. For a given value of 
the diffusion ratio, such as D = 0.5, and a fixed con- 
version rate j3 c = 0.5 the tricritical line is only slightly 
shifted within the i?-/3-plane while varying the density 
parameter a over the permitted interval (a £ [0,1]) 
as may be concluded from Fig. [SJ Traveling waves al- 
ways bifurcate supercritically beneath the tricritical line 
(dashed and dotted lines in Fig. |SJ and a change to a 
subcritical bifurcation takes place while approaching the 
codimension-2 line visualized by the solid line in Fig. [S] 
Quite a similar behavior has been found in other sys - 
tems, for instance in binary fluid convection [50l. l5ll l52| . 
If D is in- or decreased compared to D — 0.5 in Fig. [S] 




FIG. 7: The coefficient gi appearing in the amplitude 
equations Eqs. I|38fl is drawn as a function of the conver- 
sion rate ratio (5 for different diffusion ratios D, with D — 
0.3, 0.4, 0.5, 0.6, 0.7, 0.8 from left to right. The zero points il- 
lustrate how the tricitical point is moving for a given charge 
ratio 7? = —5.5 to increasing values of /3. Further parameters 
are /3 C = 0.5, a = 0.8. 
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0.26 0.32 0.38 0.44 ^ 0.50 

FIG. 8: The codimension-2 line (solid line) as well as the 
tricritical lines are plotted in the i?-/9-plane for the density 
parameters a — 0.2 (dashed line) and a = 0.8 (dotted line). 
Furthermore D = 0.5 and j3 c = 0.5. 



the tricritical lines are moved along the corresponding 
codimension-2 lines towards more negative charge ratios 
while simultaneously abuting even further to the asso- 
ciated codimension-2 lines. Thus a broadening of the 
domain of supercritically bifurcating waves is observed 
for slightly negative charge ratios. 



2. Stability of traveling waves in absence of the large-scale 
field C 

Among the classes of solutions admitted by the com- 
plex Ginzburg-Landau equation Eq. (|43|l is the family of 
plane wave solutions 




(44) 
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0.0 0.5 1.0 1.5 2.0 2.5 



FIG. 9: The positions of the codimension-2 line (solid line), 
the tricritical line (dashed line) and the line caracterizing 
marginally long-wavelength stable plane waves (dotted line) 
are displayed in the i?-/3-plane. Traveling waves bifurcate 
supercritically in the area located underneath the tricritical 
line and plane wave solutions given by Eq. 14411 are long- 
wavelength unstable in the region lying in-between the tri- 
critical line and the dotted line, wherein Co ci < —1 holds. 
Parameters used are D = /3 C = 0.5, as well as a = 0.8. 



1 + ClCQ 




0.00 0.06 0.12 0.18 0.24 



FIG. 10: The Benjamin-Feir criterion 1 + ci Co = for 
marginally long-wavelength stable plane wave solutions is 
plotted as a function of the conversion rate ratio j3 while as- 
suming a slightly negative charge ratio R — —0.8. Along the 
directed arrow the density parameter a is progressively in- 
creased from 0.2 to 0.8 by steps of 0.1 and evinces that for 
a — > 1 the domain of unstable plane waves lying in the im- 
mediate neighborhood of the codimension-2 line is incessantly 
aggrandized. Remaining parameters read D = f3 c — 0.5. 



with the frequency 

O = -a - ?yci + (cq - ci)q 2 (45) 

existing for q 2 < rj in case of a supercritical bifurca- 
tion. One easily sees that solutions of the type 144(1 
are long-wavelength stable for a range of wave numbers 
surrounding the homogeneous q = state as long as 
the Benjamin-Feir criterion, 1 + cqCi > 0, is satisfied 

0, M M El E3- 

Provided that the charge ratios R are sufficiently 
negative, traveling waves are, according to Fig. lin- 
early stable in a large area of the i?-/3-plane under- 
neath the dotted line, determined by the marginal sta- 
bility criterion 1 + c\Cq = 0, and become unstable in- 
between the latter line and the tricritical one visualized 
by the dashed line. Additionally a small unstable stripe, 
wherein the Benjamin-Feir criterion is violated, emerges 
for slightly negative charge ratios in the neighborhood of 
the codimension-2 line depicted by the solid line in Fig. [5] 
wherein D = (3 C = 0.5 and a — 0.8. 

Whereas the extension of the instability region at 
strongly negative charge ratios is barely affected by the 
choice of the density parameter a, the one located at 
mildly negative values of R utterly depends on a as 
shown exemplarily in Fig. E3f° r R = —0.8, the diffusion 
ratio D as well as the conversion rate [3 C still evaluating 
to 0.5. The marginal stability condition 1 + c\Cq = be- 
ing plotted in Fig. llOl as a function of the conversion ratio 
/? up to the associated codimension-2 point, it turns out 
that this unstable stripe vanishes for small density pa- 
rameters. If a is in contrast kept constant at 0.8 while 
the diffusion ratio D is varied, Fig. ^2p rov es that a dif- 
fusion ratio of D = 0.5 induces a maximal broadening 
of this very unstable band situated at slightly negative 



1 + CiC() 




0.00 0.08 0.16 0.24 



FIG. 11: The dependance of the marginal Benjamin-Feir cri- 
terion 1 + Ci Co = on the diffusion ratio D as a function of 
j3 for a density parameter a = 0.8, the charge ratio R eval- 
uating to —0.8. Increasing the diffusion ratio from D = 0.3 
via D = 0.4 to D — 0.5 along the arrow leads to a maximal 
enlargement of the unstable regime of wave solutions, the lat- 
ter effect being once again reversed if D > 0.5. Furthermore 
p c = 0.5. 



charge ratios. For strongly negative charge ratios the re- 
gion of Benjamin-Feir-stable traveling waves is widened 
in a similiar manner to the one of supercritically bifur- 
cating waves presented in Fig. [7| 

The effect that traveling waves become unstable via 
Benjamin-Feir resonance, which may eventually lead to 
spatiotemporal chaos, while approaching the tricritical 
point lying in the vicinity of a codimension-2 point has 
been predicted in a general context |56| and discussed 
within the scope of binary fluid convection |57| . 
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C. Stability of traveling waves in presence of the 
large-scale mode C 

In the present section we consider once more traveling 
wave solutions of the form Eq. I|44|l and study, in contrast 
to Sec. IIV B 21 their stability by implicating the long- 
wavelength stationary excitation C, that is by taking into 
account the global conservation law imposed upon the ion 
channel densities. To analyse the stability of the afore- 
said traveling waves we set 

A = [F + Vl (x,t)}e i{qx - nt) , (46a) 
C = v 2 {x,t), (46b) 
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substitute these expressions into the amplitude equations 
Eqs. (|38[) and linearizing with respect to the small per- 
turbation terms vi(x,t) and v 2 (x,t) yields in case of the 
homogeneous q = state the following coupled partial 
differential equations 

r[d t + v g d x - iO]ui = [(// + ia) + + ic a )dl] vi 

+ F [55(1 + ic 5 )d x - g 2 (l + ic 2 )] v 2 
- 5 i (1 + icx){v\ + 2 Vl )F 2 , (47a) 
d t v 2 = Dtfl v 2 - bFd x { Vl + vt) (47b) 

with F and Q denoting amplitude respectively frequency 
of the plane wave solution Eq. (|44|l . It is beneficial to 
consider sideband instability by setting 

vt(x, t) = Vi e st+iKx + V 2 * e ^ t - iKx , (48a) 
v 2 (x, t) = V 3 e ^ t+iKx + v*e^ t - iKx , (48b) 

where Vi (i = 1,2,3,4) designate complex amplitudes. 
From the coupled equations Eqs. I|47|l the aforementioned 
ansatz then generates a set of homogeneous, linear equa- 
tions whose solubility condition leads to the growth rate 
X(if) as a function of the perturbation wave number K. 

The expansion of the latter growth rate with respect 
to K 

f 2 

E(K) tt-iv g K-{l + cod) ^K 2 + (K*) (49) 

T 

shows that in the limit of small perturbation wave num- 
bers K the Benjamin-Feir criterion of Sec. IIV B~2l rcmains 
always unchanged when taking into account the large- 
scale mode C and switching to a frame moving with 
group velocity v g . Since the large-scale mode C leaves 
the stability of traveling waves untouched with respect to 
long-wavelength perturbations, it might alter their sta- 
bility by triggering an instability at finite perturbation 
wave number which is assured by Fig. Iltjf fe'l . 

Compared to the case C = studied in Sec. lIVB~2l ex- 
emplarily for D — (3 C — 0.5 and a — 0.8 through Fig. [5J 
Fig.^]clearly evinces a strong influence of the large-scale 
mode C due to the conservation law on the stability of 
traveling waves: Regions of the i?-/3-plane, wherein sta- 
ble plane wave solutions occur, are significantly restricted 



FIG. 12: Stability regions of nonlinear traveling waves given 
by Eqs. lUSl within the i?-/3-plane for D = f3 c = 0.5 and 
different values of the density parameter a. The solid line vi- 
sualizes the codimension-2 line, traveling waves being stable 
in the area lying underneath the lines of marginal stability, 
namely the dashed-dotted one corresponding to a = 0.6, the 
dashed one to a = 0.75 and the dotted one to a density pa- 
rameter of a = 0.8. The distance to the bifurcation point was 
chosen to be r\ — 0.03. 

and shifted towards more negative charge ratios. Fig. 1121 
shows those stability domains limited, depending on the 
choice of the density parameter a, to the area located 
below the dashed-dotted, the dashed or the dotted line 
of marginally stable nonlinear traveling waves satisfying 
£(if) = 0. The real part of the latter growth rate T,(K) 
is plotted in Fig. I13f b) as a function of the perturba- 
tion wave number K presuming the density parameter a 
evaluates to 0.8. In account with the stability diagram, 
i.e. Fig.[T21 the tuple (i? = -5.5, /3 = 2.17) corresponds 
to wave solutions becoming unstable with respect to fi- 
nite perturbation wave numbers K as expected from the 
power series of £(K) performed in Eq. (|49|) . To hight- 
light the effect of the large-scale mode C evoked by the 
conservation law, the real part Re(S) of the same growth 
rate is depicted in Fig. EH a ) hi case °f a vanishing real 
mode C. The stability analysis of C = in Sec. IIV B 21 
has shown by means of the stability diagram Fig. [S] that 
for the parameter set (R = —4, /3 = 1.86, a = 0.8) trav- 
eling waves are long-wavelength unstable and thus the 
observed instability in Fig. I13f a) is generically different 
from the one encountered in the presence of a conserved 
quantity whose typical shape is given by Fig ll3f bh 

Moreover the sensitivity of the former unstable stripe 
at slightly negative values of R (cf. enlargement in Fig.^J) 
on the density parameter, depicted in Fig. 1101 is ampli- 
fied as it may be concluded from Fig. ^] as well as from 
Fig. 1141 wherein the points of marginally stable long- 
wavelength nonlinear plane wave solutions are plotted as 
a function of a for a constant conversion ratio (3. The lat- 
ter figure also stipulates that, for a given value of /3, there 
exists a certain density parameter a generating a maxi- 
mal outspread of the area of stable traveling waves. For 
the given parameter set D = (3 C = (3 — 0.5 and rj = 0.03 
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FIG. 13: The typical shape of the real part Re(E) corre- 
sponding to the growth rate T,(K) specified by the ansatz 
Eq. I|48fl is plotted as a function of the perturbation wave 
number K for D = /3 C = 0.5, a — 0.8 as well as r\ — 0.03 in 
(b). R = —5.5 and (3 = 2.17 leads, according to the stability 
diagram in Fig. 1121 to wave solutions unstable with respect 
to finite perturbation wave numbers K. On the left hand 
side however the tuple (R = — 4, f3 = 1.86) induces, in the 
absence of the large-scale field C characterizing the conserva- 
tion law, generic long-wavelength unstable traveling waves in 
accordance with Fig.|2l 




FIG. 14: Points of marginal stability for a fixed conversion 
ratio f3 = 0.5 as a function of the density parameter a, the 
diffusion ratio D and the conversion rate /3 C evaluating to 0.5 
as in Fig. El Furthermore i) = 0.03. 



a density parameter of a ~ 0.65, meaning the ion chan- 
nels contribute up to 65% to the total conductance of 
the membrane layer, realizes the largest extension of the 
stable traveling wave regime. 



D. Numerical results 

A major problem always encountered when dealing 
with amplitude equations is that their validity range 
around the threshold r\ = remains a priori unpre- 
dictable. Therefore we have determined the amplitude 
of a stripe solution by numerical simulations of the basic 
Eqs. (|llfl as a function of the relative distance r\ to the 
bifurcation point and compared these numerical results 
with the analytical solution A = \Jr\j g\ eo (N ,N C ) T in 
Fig. Up to r\ ~ 0.05 we find for the particle density 
N c of the closed ion channels a fairly good agreement be- 



FIG. 15: The amplitudes N and N c of traveling (opened 
and closed respectively) ion channel densities obtained by nu- 
merical simulations (crosses) of the microscopic model given 
by Eqs. 11 U are compared to the analytical solution A — 
\J r\jg\ eo (iV , N C ) T . The parameters used are D = l,/3 = 
0.4, R = —1.8, a = 0.8, t v = 0.01 with the individual conver- 
sion rates being (3 C = 0.5 and f3 = 0.2. 



tween analytical and numerical results whereas the same 
holds up to -q ~ 0.03 for the opened ion channel density 
N a , which gives a reasonable estimate of the quantitative 
validity range of the preceding results obtained in terms 
of amplitude equations. 



V. DISCUSSION AND CONCLUSION 

Besides its biological motivation, the model presented 
in this work exhibits interesting generic properties con- 
cerning pattern formation that are for the first time 
discussed within the scope of a specific model system. 
The homogeneous ion channel density distribution on the 
biomembrane may become unstable with respect to ei- 
ther a stationary or an oscillatory bifurcation, whereas 
we focused exclusively on analyzing the latter bifurca- 
tion type. As the total ion channel density is a conserved 
field, the supercritical Hopf-bifurcation occuring at finite 
wave number shows nonlinear properties largely differing 
from those emerging in systems lacking conserved order 
parameters. 

In presence of a global conservation law, such as the 
one encountered in our model, it can be expected that the 
asymptotic description of pattern formation by means of 
amplitude equations is altered through the coupling of 
known forms of envelope equations to long-wavelength 
modes. This allegation has been confirmed by the per- 
turbation expansion of the initial model equations that 
yielded a set of coupled amplitude equations in case of 
traveling ion channel density waves rather than a single 
complex Ginzburg-Landau equation. The validity range 
of the latter reduced dynamics has been checked by sim- 
ulations of the inital model equations. This coupling to 
stationary, long- wavelength excitations has however been 
completely neglected during the discussion of the station- 
ary bifurcation scenario in Ref. |43| which gives conse- 
quently a distorted picture of the nonlinear bifurcation 
behavior. 

Detailed phase diagrams highlight in which param- 
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eter space Hopf-bifurcations at finite wavelengths take 
place. Neither the threshold value of the control param- 
eter nor the supercriticality of the oscillatory bifurcation 
are modified by implicating the long-wavelength mode 
evoked by the conserved ion channel density field into 
the asymptotic description near onset of instability. The 
stability of traveling ion channel density waves is however 
strongly affected by the latter real mode for infinitely ex- 
tended systems. It turns out that the long-wavelength 
mode, reflecting the conserved density field, triggers an 
instability at finite perturbation wave number whereas 
plane wave solutions would be long-wavelength unstable 
in the absence of the large-scale mode. If compared to a 
case without a real mode, stable traveling wave solutions 
are generally shifted towards strongly negative ratios of 
opened and closed ion channel charges respectively and 
thus the domain, wherein spatiotemporal chaos is likely 
to be seen, significantly broadens. If finally all model pa- 
rameters are fixed the observation of stable traveling ion 
channel density waves becomes less probable the more 
the ion channels contribute to the total conductance of 
the considered membrane layer. 

The dynamics of traveling waves in finite geome- 
tries and consequently their stability properties natu- 
rally differ from those in infinitely extended systems 
[HsL IHol IfioL l6ll ] . How the finite size of the system fur- 
ther influences the studied traveling ion channel density 
waves at hand, already affected by a large-scale mode, is 
postponed to future work. 

The formation of dissipative patterns of ion channels 
on a membrane that were discussed in the present work is 
most likely to be observed in in vitro systems as the one 
proposed in Ref. where a linear array of 48 field-effect 
transistors in silicon was coated with a bimolecular layer 
of lipid and gramicidin. Nevertheless the investigated 
electrodiffusive mechanism seems to be also relevant for 
in vivo systems as shown by experiments on the effect of 
electric fields on clustering of acetylcholine receptors as 
cited in Ref. H|. 

Couplings between a long- wavelength mode and a trav- 
eling wave arc expected in other models from cell biology 
as well, e.g. in mixtures of cytoskeletal filaments inter- 
acting with molecular motors [64[. The investigation of 
an oscillatory instability occuring at finite wave num- 
ber coupling to a stationary, long-wavelength mode, that 
was motivated by the biological model at hand, will be 
extended within the scope of a reaction-diffusion system 
in a forthcoming work to oscillatory, long-wavelength ex- 
citations influencing on oscillatory instabilities at finite 
wave number 491. 
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APPENDIX A: AMPLITUDE EXPANSION 

According to the slow time and spatial scales T\, T 
and X introduced in Sec. II V Al all time and space deriva- 
tives occuring in the scaled equations Eqs. (|TTT) may be 
substituted by 

dt~> d t + y/rj d Tl +r)d T , (Ala) 
d x ^d x + ^nd x . (Alb) 

The pattern-generating ion channel densities N a and N c 
as well as the transmembrane voltage V have to be ex- 
panded with respect to powers of y/rj 

No^VvNl+v N 2 + 17* N 3 + O (r/ 2 ) , (A2a) 
N c = y/rj Nl + r) N 2 + v i N 3 + O (v 2 ) , (A2b) 
V= y/rjV 1 +r}V 2 + r}?V 3 + 0(r} 2 ) (A2c) 

in order to accomplish a perturbation expansion of 
Eqs. JUJ. Within the adiabatic limit, i.e. if instanta- 
neously adapting membrane voltage V is assumed, mean- 
ing formally C m = or t v — 0, Eq. Ijllcfl becomes time 
independant and may be solved successively by means of 
the resulting equation hierarchy in V 1 (i = 1, 2, 3) 



: TV 1 = 
775 : TV 2 = 
r)i : TV 3 = 



-e c Nl , (A3a) 
-e c 7V 2 + (2d x 0x - aNDV 1 , (A3b) 
-e^ + N 3 ) + (d 2 x - aN^V 1 
F (2d x d x ~ aN l )V 2 , (A3c) 



yielding the three leading order contributions to V 



V 1 
V 2 
V 3 



F- 1 [~e»Nl] 



H N 2 + (2d x dx-aN 1 )V 1 ] 



e 

(2d x d x 



O 2 



X 



aN^V 1 



aNl)V 2 



wherein the operator 



T = 1 - di 



(A4a) 
(A4b) 

(A4c) 
(A5) 



has been introduced. 

Since the desired solutions iV* and N l c with i = 1, 2, 3 are 
harmonic functions of e *("ct±fcca;) ; ^n e occuring operator 
J~ in the antecedent expressions for V 1 may be replaced 
in the forthcoming calculations as 



a ni(k c x^-uj c t) 



1 



1 + n 2 k 2 



ni(k c x-\-uj c t) 



(A6) 



with n = 0,1,2, ... 

By performing the perturbation expansion of the 
scaled equations Eqs. Qllafl . Ijllb|) governing the spa- 
tiotemporal evolution of the ion channel densities and 
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taking into account Eqs. (|A4|) as well as Eq. I|A6(1 one wherein the occuring differential operators read 
ends up with a hierarchy in ^Jr\ 

N °i), (* = 1,2,3), (A8) 

rp :£ w 2 =M(wi)-(A+>to&r 1 )wi (A7b) /l0\ 

3 Mo =Mi = M 2 = r, i , (A9) 

?7 2 :£ow 3 = .A/ 2 (wi,w 2 ) - (A w 2 V u 1 / 

- (£ 2 + M 2 <9 T ) wi (A7c) 



4> = 



£ 2 = 



Ci = 2 



AA 2 



-d\ + ef + fl* ) ^ 2 + 53% ^ + 4d%] ^ 3 
DQe? [(dl + d 2 x ) T 2 + 5d 2 d 2 x T + 4d%] -Dd 

( 

-T-'dl{NlT-'Nl)]ae^- 




d x d x , 



-D 




d x NlT^Nl+NlT-^d 2 x Nl 
\ DQ [d x Nl F-i-Nl + Nl T-'dlNl] 
T- X dl {N 2 T^Nl) - 2T- 1 d x d x (l + d 2 ) (Nl T-'Nl) - mi] 



m 2 d x Nl + (d x Nl + d x N 2 + N 2 d x ) T^d^l + m 3 N l 
DQ [m 2 d x Nl + (d x Nl + 8 X N 2 + N 2 d x ) T^d^l + m 3 Nl] 



with the abbreviations 



mi = - d x T- 1 [Nl {T- X N 2 + 2d x d x T- 2 Nl - a^ 1 (JV&^iV*))] , 

m 2 = (2d x d x T- 2 + dxT- 1 ) Nl + T^N 2 , 

m 3 = (2d 3 x d x T- 2 + 2d x d x T- 1 ) N X Q + d^N 2 . 

I 



The aforementioned hierarchy, Eqs. l|A7j) . has to be 
solved successively for Wj with i = 1,2,3. The first of 
these three equations is fulfilled if the ansatz given in 
Eq. i|35|l is chosen for wi. The determined solution may 
be inserted into equation Eq. l|A7b(l which leads to a lin- 
ear, inhomogeneous equation of the type £o w 2 = 1 2 - 
Since Cq is not invertible, as it has a zero eigenvalue, the 
only possibility for the latter equation to have nontriv- 
ial solutions is that I 2 is orthogonal to the kernel of £q 
which yields the solvability condition 

(foe^-r+^l^) (All) 
= i jfdrl I_dte-^^4X 2 = 



a second Fredholm's alternative 



(z |Z 2 > =0 



(AlOa) 



(AlOb) 



(AlOc) 



(AlOd) 



(AlOe) 
(AlOf) 
(AlOg) 



(A12) 



where zq = (1, —(3) designates the left eigenvector span- 
ning the adjungated kernel of £o\ k=0 - Inserting the so- 
lution w 2 of Eq. (|A7b|> as well as Wi into Eq. (|A7c|) 
gives rise to another linear, inhomogeneous equation of 
the form £o w 3 — ^3- To the latter equation are associ- 
ated two solvability equations 



(f oe i(k?- r +-ct)| X3) = 0) 

(zol^s) =0. 



(A13) 
(A14) 



with f denoting the left eigenvector of £o ■ Terms propor- 
tional to e°, linked to the large-scale mode C, stipulate 



From Fredholm's alternatives at order 77, namley 
Eqs. (|Alip and HA12(1 . as well as at order ^/ry 3 , that is 
Eqs. (|A13|I and IjAMfl . the final amplitude equations |(38J| 
can be obtained after rescaling space and time according 
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to as suggested by IjAll) . 

d t A = ^d Tl A + v d T A, (A15) 
d t C = y/Tjd Tl C + r)d T C (A16) 
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